function [H,D,Psum] = power_distribution(G)
%power_distribution 此处显示有关此函数的摘要
%   G——一个有向图，每条边的方向代表有功功率流的方向。
%       Edges应具有如下字段：Endnodes（图自带的），P1（首端送出功率，正数），P2（末端接收功率，正数）
%       Nodes应具有如下字段：Pg（每个节点发电机组/上级电网的功率，正数），Pd（每个节点负荷功率，正数）
%   H——节点输出分布因子矩阵
%   D——发电机至节点路径输出分布因子矩阵
%   Psum——有功通量
%   
Nn=numnodes(G);
Ng=sum(G.Nodes.Pg>0);
Gloc=find(G.Nodes.Pg>0);
Ne=numedges(G);

%% 计算H
%有功通量部分
Psum=zeros(Nn,1);
idx=G.Nodes.Pg>0;
Psum(idx)=G.Nodes.Pg(idx);

for k=1:Ne
    s=G.Edges.EndNodes(k,[1 2]);
    Psum(s(2))=Psum(s(2))+G.Edges.P2(k);
end

H.ii=zeros(1,Ne);
H.jj=zeros(1,Ne);
H.val=zeros(1,Ne);
for k=1:Ne
    s=G.Edges.EndNodes(k,[1 2]);
    H.ii(k)=s(1);
    H.jj(k)=s(2);
    H.val(k)=G.Edges.P2(k)/Psum(s(1));
end

%对角线部分
H.ii=[1:Nn H.ii];
H.jj=[1:Nn H.jj];
H.val=[(G.Nodes.Pd./Psum).' H.val];

%转换为稀疏矩阵
H=sparse(H.ii,H.jj,H.val);

%% 计算D
A=sparse(diag(2*ones(1,14)))-H.';
D=zeros(Ng,Nn);
tol=1e-10;
for k=1:Ng
    Ag=A;
    Ag(Gloc(k),Gloc(k))=Ag(Gloc(k),Gloc(k))+1;
    Ag(Gloc(k),:)=Ag(Gloc(k),:)-A(Gloc(k),:);
    b=zeros(Nn,1);
    b(Gloc(k))=1;
    [x,flag,relres] =lsqr(Ag,b);
    x(x>-tol & x<tol)=0;
    D(k,:)=round(x.',5);
end

